{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d92ee154",
   "metadata": {},
   "source": [
    "# Readout Cavity Calibration\n",
    "*Copyright (c) 2021 Institute for Quantum Computing, Baidu Inc. All Rights Reserved.*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "230b6f9b",
   "metadata": {},
   "source": [
    "## Outline\n",
    "This tutorial introduces the simulation of readout cavity calibration using the readout simulator. The outline of this tutorial is as follows:\n",
    "- Introduction\n",
    "- Preparation\n",
    "- Calibrating the Readout Cavity Transition Frequencies\n",
    "- Calibrating the Dispersive Shift and Coupling Strength\n",
    "- Measuring the decay rate\n",
    "- Summary"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b6b5b3e",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "In superconducting circuit, to acquire the state of a qubit, we can probe the readout cavity coupled with this qubit to achieve the qubit state indirectly. Concretely, we first apply readout pulse signal and then detect and analyze the reflected signal. Because the phase shift and amplitude change depend on the qubit state, we are able to know whether the outcome is \"0\" or \"1\" by this change.\n",
    "\n",
    "In the real experiment of calibration, the first step to is to find the parameters of readout cavity. This tutorial introduces how to use Quanlse to simulate the readout cavity calibration.\n",
    "\n",
    "A coupled cavity-qubit system can be described by the Jaynes-Cummings Hamiltonian in the dispersive regime \\[1\\]:\n",
    "\n",
    "$$\n",
    "\\hat{H}_{\\rm JC} = \\omega_r \\hat{a}^\\dagger \\hat{a} + \\frac{1}{2}\\omega_q \\hat{\\sigma}_z + \\chi \\hat{a}^\\dagger \\hat{a} \\hat{\\sigma}_z,\n",
    "$$\n",
    "\n",
    "\n",
    "where $\\hat{a}$, $\\hat{a}^\\dagger$ are annihilation and creation operators and $\\hat{\\sigma}_z$ is the Pauli-Z operator. $\\omega_r$ and $\\omega_q$ denote the bare frequencies of the readout cavity and the qubit, $\\chi$ is the dispersive shift and takes the form \\[2\\]:\n",
    "\n",
    "$$\n",
    "\\chi = \\frac{g^2 \\alpha}{\\Delta_{qr}(\\Delta_{qr} + \\alpha)}.\n",
    "$$\n",
    "\n",
    "where $\\alpha$ is the qubit anharmonicity, $\\Delta_{qr} = \\omega_q - \\omega_r$ is qubit-cavity detuning and $g$ is the qubit-cavity coupling strength. The interaction term $\\chi \\hat{a}^\\dagger \\hat{a} \\hat{\\sigma}_z$ in $\\hat{H}_{\\rm JC}$ gives rise to a shift of $2\\chi$ in the transition frequency of the readout cavity when the qubit state is $|0\\rangle$ and $|1\\rangle$. Therefore in the experiment, by performing frequency sweep for the cavity with qubit state prepared in $|0\\rangle$ or $|1\\rangle$ respectively, we can obtain transition frequency $f_0$ and $f_1$ and therefore the frequency shift $2\\chi$. Finally, the cavity-qubit couping strength is indirectly calculated using the expression above. \n",
    "\n",
    "In addition to the transition frequency and dispersive shift, the linewidth $\\kappa$ can be measured to determine the photon decay rate of the readout cavity. To simulate the interaction between cavity-qubit system and the environment, the evolution of the system density matrix $\\hat{\\rho}(t)$ is given by Lindblad master equation \\[3, 4\\]:\n",
    "\n",
    "$$\n",
    "\\frac{d \\hat{\\rho}(t)}{dt} = -i[\\hat{H}(t), \\hat{\\rho}(t)] + \\frac{\\kappa}{2}[2 \\hat{a} \\hat{\\rho}(t) \\hat{a}^\\dagger - \\hat{\\rho}(t) \\hat{a}^\\dagger \\hat{a} - \\hat{a}^\\dagger \\hat{a} \\hat{\\rho}(t)].\n",
    "$$\n",
    "\n",
    "\n",
    "The decay rate is therefore acquired by fitting the spectrum and extract the linewidth in the Lorentzian function.\n",
    "\n",
    "The observable quantities we take are the two orthogonal quadratures $\\hat{X} = \\frac{1}{2}(\\hat{a}^\\dagger + \\hat{a})$ and $\\hat{Y} = \\frac{i}{2}(\\hat{a}^\\dagger - \\hat{a})$. In the experiment, through a series of signal processing on the pulse reflected from the readout cavity, we can obtain the voltage $V_I$ and $V_Q$ related to these two orthogonal quadratures.\n",
    "\n",
    "In this tutorial, we simulate the calibration of readout cavity by solving the qubit-cavity dynamics: the cavity transition frequencies with different qubit states ($|0\\rangle$ and $|1\\rangle$) $\\omega_{r0}$ and $\\omega_{r1}$, the linewidth $\\kappa$ and the dispersive shift $\\chi$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b7b6f28f",
   "metadata": {},
   "source": [
    "## Preparation\n",
    "\n",
    "To run this tutorial we need to import the following necessary packages from Quanlse and other commonly-used Python libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "774fbc17",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import tools from Quanlse\n",
    "from Quanlse.Superconduct.Simulator.ReadoutSim3Q import readoutSim3Q\n",
    "from Quanlse.Superconduct.Calibration.Readout import resonatorSpec, fitLorentzian, lorentzian\n",
    "\n",
    "# Import tools from other python libraries\n",
    "from scipy.signal import find_peaks\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from math import pi"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "faeffc92",
   "metadata": {},
   "source": [
    "## Calibrating the Readout Cavity Transition Frequencies\n",
    "\n",
    "In this section, we will calibrate the transition frequencies of the reading cavity when the qubit is in the ground state and the first excited state, respectively. Here, we first use the predefined function `readoutSim3Q()` to return the Class object `readoutModel` containing information of the readout cavity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3cff0622",
   "metadata": {},
   "outputs": [],
   "source": [
    "readoutModel = readoutSim3Q()  # Initialize a readoutModel object"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9cca0cd",
   "metadata": {},
   "source": [
    "Then, we set the range of frequency sweep `freqRange`, the drive amplitude `amp` and the duration of the readout pulse `duration`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f232330f",
   "metadata": {},
   "outputs": [],
   "source": [
    "freqRange = np.linspace(7.105, 7.125, 300) * 2 * pi  # the range of frequency to probe the resonator, in 2 pi GHz\n",
    "amp = 0.0005 * (2 * pi)  # drive amplitude, in 2 pi GHz\n",
    "duration = 1000  #  duration of the readout pulse, in nanoseconds"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7b075c8",
   "metadata": {},
   "source": [
    "Use the function `resonatorSpec` to simulate the frequency sweep of the readout cavity when qubit is in the ground state, and input the index of the resonator `onRes`, the range of the frequency sweep `freqRange`, the amplitude `amp` and the duration `duration` with `qubitState` set to be in the ground state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1fc95a43",
   "metadata": {},
   "outputs": [],
   "source": [
    "vi0, vq0 = resonatorSpec(readoutModel=readoutModel, onRes=[0], freqRange=freqRange, \n",
    "                         amplitude=amp, duration=duration, qubitState='ground')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c3591f5a",
   "metadata": {},
   "source": [
    "The result returns the measured signal $V_I$ and $V_Q$. We plot $V_Q$ (or $V_I$) with respect to the drive frequency."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5612aa6e",
   "metadata": {},
   "outputs": [],
   "source": [
    "idx0 = find_peaks(vq0[0], height=max(vq0[0]))[0]  # find the index of the transition frequency\n",
    "w0 = freqRange[idx0][0]  # transition frequency\n",
    "print(f'The resonator transition frequency with qubit in ground state is {(w0 / (2 * pi)).round(3)} GHz')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6f454dee",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(freqRange / (2 * pi), np.array(vq0[0]))\n",
    "plt.plot()\n",
    "plt.xlabel('$\\omega_d$ (GHz)')\n",
    "plt.ylabel('signal (a.u.)')\n",
    "plt.title('Readout resonator spectrum')\n",
    "plt.vlines((freqRange / (2 * pi))[idx0], 0, max(vq0[0]), linestyles='dashed')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b03e96f",
   "metadata": {},
   "source": [
    "From the result of the simulation shown above, we can see that the read cavity transition frequency is around 7.118 GHz when the qubit is in the ground state. Next, we calibrate the read cavity transition frequency when the qubit is in the excited state using the same procedure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23f8e3f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "vi1, vq1 = resonatorSpec(readoutModel=readoutModel, onRes=[0], freqRange=freqRange, \n",
    "                         amplitude=amp, duration=duration, qubitState='excited')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2c4f7fa7",
   "metadata": {},
   "outputs": [],
   "source": [
    "idx1 = find_peaks(vq1[0], height=max(vq1[0]))[0]\n",
    "w1 = freqRange[idx1][0]\n",
    "print(f'The resonator transition frequency with qubit in excited state is {(w1 / (2 * pi)).round(3)} GHz')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a99db7fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(freqRange / (2 * pi), np.array(vq1[0]))\n",
    "plt.plot()\n",
    "plt.xlabel('$\\omega_d$ (GHz)')\n",
    "plt.ylabel('signal (a.u.)')\n",
    "plt.title('Readout resonator spectrum')\n",
    "plt.vlines((freqRange / (2 * pi))[idx1], 0, 2, linestyles='dashed')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb3c9c38",
   "metadata": {},
   "source": [
    "It can be seen in the spectrum that the readout cavity transition frequency is about 7.112 GHz when the qubit is in the first excited state."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "642c794e",
   "metadata": {},
   "source": [
    "## Calibrating the Dispersive Shift and Coupling Strength \n",
    "\n",
    "In the previous section, we obtained the calibrated frequencies $f_0$ and $f_1$, so that the dispersion shift $\\chi$ can be calculated directly by,\n",
    "\n",
    "$$\n",
    "\\chi = \\frac{|f_0 - f_1|}{2}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "182be2e8",
   "metadata": {},
   "outputs": [],
   "source": [
    "chi = abs(w0 - w1) / 2\n",
    "print(f'The dispersive shift is {(chi * 1e3 / (2 * pi)).round(3)} MHz')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2deb56e6",
   "metadata": {},
   "source": [
    "Combining the expressions of $\\chi$ given in the \"Introduction\" section, we can derive the expression of cavity-qubit coupling strength in terms of other known parameters:\n",
    "\n",
    "$$\n",
    "g = \\sqrt{\\frac{\\chi\\Delta_{qr}(\\Delta_{qr}+\\alpha)}{\\alpha}}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3dd27411",
   "metadata": {},
   "source": [
    "Extract the theoretical parameters from `readoutModel` and calculate the coupling strength $g$ given above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "67e77118",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract parameters from the model\n",
    "g = readoutModel.coupling[0]  # therotical qubit-resonator coupling strength \n",
    "wq = readoutModel.pulseModel.qubitFreq[0]  # qubit bare frequency\n",
    "alpha = readoutModel.pulseModel.qubitAnharm[0]  # qubit anharmonicity\n",
    "wr = (w0 + w1) / 2  # estimated resonator frequency\n",
    "detuning = wq - wr  # qubit-resonator detuning"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e816795b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# coupling strength calculation\n",
    "\n",
    "def qrCoupling(chi, detuning, alpha):\n",
    "    g = np.sqrt(abs(chi * detuning * (detuning + alpha) / alpha))\n",
    "    return g"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "718b16aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "gEst = qrCoupling(chi, detuning, alpha)  # Estimated qubit-resonator coupling strength"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "885317d0",
   "metadata": {},
   "source": [
    "Compare the theoretical value and the estimated value of $g$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "32638ebb",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f'Theoretical coupling strength is {g * 1e3 / (2 * pi)} MHz')\n",
    "print(f'Estimated coupling strength is {(gEst * 1e3 / (2 * pi)).round(1)} MHz')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b5abb80e",
   "metadata": {},
   "source": [
    "The coupling strength of the readout cavity and the qubit, obtained by calibrating the dispersion shift and indirect calculations, is 132.4 MHz, which is in good agreement with the theoretical value of 134.0 MHz."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd0e915e",
   "metadata": {},
   "source": [
    "## Measuring the decay rate\n",
    "\n",
    "After we have the spectrum of cavity frequency, we are able to estimate the decay rate $\\kappa$ by the linewidth of the Lorentzian function. Here, we use the function `fitLorentzian`, input the range of frequency sweep and the reflected signal to fit the spectrum and estimate the linewidth $\\kappa$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "91e8b986",
   "metadata": {},
   "outputs": [],
   "source": [
    "param, cov = fitLorentzian(freqRange, vq0[0])  # Fit the curve using lorentzian function\n",
    "kappaEst = param[2]  # Estimated linewidth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "98d2f1f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(freqRange / (2 * pi), lorentzian(freqRange, param[0], param[1], param[2], param[3]), '.')\n",
    "plt.plot(freqRange / (2 * pi), vq0[0])\n",
    "plt.xlabel('$\\omega_d$ (GHz)')\n",
    "plt.ylabel('signal (a.u.)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "97df7010",
   "metadata": {},
   "source": [
    "Compare the theoretical value and the estimated value of decay rate (or linewidth)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6188135b",
   "metadata": {},
   "outputs": [],
   "source": [
    "kappa = readoutModel.dissipation\n",
    "\n",
    "print(f'Theoretical linewidth is {kappa * 1e3 / (2 * pi)} MHz')\n",
    "print(f'Estimated linewidth is {(kappaEst * 1e3 / (2 * pi)).round(3)} MHz')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "210d12e9",
   "metadata": {},
   "source": [
    "From the simulation results, we can see that the decay rate $\\kappa$ set in the master equation is 2.0 MHz, while the linewidth obtained from the spectrum is 1.987 MHz, indicating that the interaciton strength between the reading cavity and the environment can be indirectly calibrated by scanning the frequency of the reading cavity and calculating the line width in the experiment."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2866975",
   "metadata": {},
   "source": [
    "## Summary\n",
    "Users can click on this link [tutorial-readout-cavity-calibration.ipynb](https://github.com/baidu/Quanlse/blob/main/Tutorial/EN/tutorial-readout-cavity-calibration.ipynb) to jump to the corresponding GitHub page for this Jupyter Notebook documentation and run this tutorial. You can try different hardware parameters of the readout cavity and run the codes in this tutorial to simulate the cavity calibration in the superconducting quantum computing experiment."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27aa0b72",
   "metadata": {},
   "source": [
    "## Reference\n",
    "\n",
    "\\[1\\] [Blais, Alexandre, et al. \"Cavity quantum electrodynamics for superconducting electrical circuits: An architecture for quantum computation.\" *Physical Review A* 69.6 (2004): 062320.](https://journals.aps.org/pra/abstract/10.1103/PhysRevA.69.062320)\n",
    "\n",
    "\\[2\\] [Koch, Jens, et al. \"Charge-insensitive qubit design derived from the Cooper pair box.\" *Physical Review A* 76.4 (2007): 042319.](https://journals.aps.org/pra/abstract/10.1103/PhysRevA.76.042319)\n",
    "\n",
    "\\[3\\] [Lindblad, Goran. \"On the generators of quantum dynamical semigroups.\" *Communications in Mathematical Physics* 48.2 (1976): 119-130.](https://link.springer.com/article/10.1007/bf01608499)\n",
    "\n",
    "\\[4\\] [Bianchetti, R., et al. \"Dynamics of dispersive single-qubit readout in circuit quantum electrodynamics.\" *Physical Review A* 80.4 (2009): 043840.](https://journals.aps.org/pra/abstract/10.1103/PhysRevA.80.043840)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
